Homicide Data Map Visualizations

Part 1: Map Visualizations for Cook County Homicides, 1940-1965

Loading packages

The following R packages will be needed:

library(sf)
library(dplyr)
library(tmap)
library(leaflet)
library(leafgl)
library(crosstalk)
library(sfdep)
library(ggplot2)

Loading Data

The homicide and geographic boundary data are loaded in. Note that 83 homicide data points had no geography associated with them and thus had to be dropped for spatial analysis. An additional 6 data points had no date, and had to be dropped as well.

# Load in homicide data and remove entries with missing geometry (n=83) and missing date (n=6)
homicides = read_sf("homicides_1940-1965_geocoded_v3/homicides_1940-1965_geocoded_v3.shp") %>%
  filter(!st_is_empty(.)) %>%
  filter(date != "NaT")

city_bounds = read_sf("Chicago_City_Limits-shp/Chicago_City_Limits.shp") %>%
  st_transform(st_crs(homicides))

hexagons = read_sf("hexagons/hexagons.shp") %>%
  st_transform(st_crs(homicides))

CAs = read_sf("Boundaries - Community Areas (current)/geo_export_e916a3b4-0ea3-4ca2-944a-cb060aacf705.shp") %>%
  st_transform(st_crs(homicides))

Simple Placement Map

A small amount of “jitter” is applied to points in order to make overlapping points both visible:

homicides_jitter = st_jitter(homicides, factor = 0.0001)

An interactive map of homicides:

leaflet(homicides_jitter) %>%
  addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
  addPolygons(data = city_bounds, opacity = 0.2) %>%
  addGlPoints(homicides_jitter, radius = 2, fillColor = "black")

Preliminary findings: Homicides 1940-1964 were concentrated primarily on the South and West Sides, with a smaller cluster appearing along the lake shore to the north as well. The Northwest and Southwest sides had a far lower concentration of homicides, though this may be due to lower population density at the time. On a more local scale, homicides tend to cluster around major thoroughfares, like Madison St on the West Side.

Choropleth Maps

Hexagons

Spatial count of homicides per hexagon:

# Assign each point the grid ID of its surrounding hexagon
homicides_hex_joined = st_join(homicides, hexagons, join = st_within)

# Create table of frequency by hexagon
hex_freq_table = as.data.frame(table(homicides_hex_joined$GRID_ID)) %>%
  rename(GRID_ID = Var1, hom_ct = Freq)

# Merge frequency table back into hexes
hexagons_hom_ct = hexagons %>%
  left_join(hex_freq_table, by = "GRID_ID") %>%
  mutate(hom_ct = ifelse(is.na(hom_ct), 0, hom_ct))

Viewing the hexagons:

tmap_mode("view")
tm_shape(hexagons_hom_ct) + tm_polygons(col = "hom_ct", breaks = c(0,1,5,10,50,100,Inf), title="Homicides", id="hom_ct", popup.vars = c("Grid ID:" = "GRID_ID", "Homicides:" = "hom_ct"))

Note that a custom classification scheme is used because of the extreme skew of the data. The 0 to 1 bin shows hexes where no homicide incidents occurred at all.

Community Areas

Spatial count of homicides per community area:

# Assign each point the grid ID of its surrounding CA
homicides_CA_joined = st_join(homicides, CAs, join = st_within)

# Create table of frequency by CA
CA_freq_table = as.data.frame(table(homicides_CA_joined$area_num_1)) %>%
  rename(area_num_1 = Var1, hom_ct = Freq)

# Merge frequency table back into CAs
CAs_hom_ct = CAs %>%
  left_join(CA_freq_table, by = "area_num_1") %>%
  mutate(hom_ct = ifelse(is.na(hom_ct), 0, hom_ct))

Viewing the community areas:

tmap_mode("view")
tm_shape(CAs_hom_ct) + tm_polygons(col = "hom_ct", breaks = c(0,1,5,10,50,100,500,Inf), title="Homicides", id="hom_ct", popup.vars = c("Area Name:" = "community", "Homicides:" = "hom_ct"))

The same classification scheme is used for consistency, though it should be noted that community areas are generally larger than hexagons. Thus, homicides exceed 500 in certain areas, and so an extra bin, “500 or more”, is added here.

Preliminary findings: As seen above, homicides tend to cluster in the South and West Sides. However, while placement maps are generally good at highlighting areas of high concentration, the two choropleth maps above elucidate areas of low concentration as well, such as the community area of North Park, which had 0 homicides during the entire time period. These results should be interpreted with caution, however, as they are not population standardized, and may simply be apparent due to a lack of population density.

Time Series Placement Map

First, homicide dates must be converted into date format. Then SharedData from crosstalk is used to link data between the slider and the map below:

homicides_jitter$start = as.Date(homicides_jitter$date)
shared_data = SharedData$new(homicides_jitter)

The map:

filter_slider("date", "Date", shared_data, ~start, width = "100%")

Date

leaflet(shared_data) %>%
  addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
  addPolygons(data = city_bounds, opacity = 0.2) %>%
  addCircleMarkers(radius = 1, weight = 0.2, color = "black")

Preliminary findings: The interactive map above illustrates the outward spread of homicides in Chicago over time, likely paralleled by the city’s outward growth. However, the change in concentration of homicides in the West Side is quite notable; while in 1940 homicides appear to be no more concentrated on the West Side than any other area of the city, by 1964 they are just as concentrated in that area, if not more, than on the South Side.

Time Series Choropleth Maps

Hexagons

Spatial count of homicides per hexagon for each year:

# Create table of frequency by hexagon and year
hex_freq_table_yr = as.data.frame(table(homicides_hex_joined$GRID_ID, homicides_hex_joined$year)) %>%
 rename(GRID_ID = Var1, year = Var2, hom_ct = Freq)

# Convert factors to integers/strings
hex_freq_table_yr$GRID_ID = as.character(paste(hex_freq_table_yr$GRID_ID))
hex_freq_table_yr$year = as.integer(paste(hex_freq_table_yr$year))

# Create an expanded grid of hexagons that includes each hex/year combo
hexagons_yr = expand.grid(GRID_ID = hexagons$GRID_ID, year = 1940:1964) %>%
 left_join(hexagons, by = "GRID_ID") %>%
 st_as_sf()

# Merge frequency table back into hex/year grid
hexagons_hom_ct_yr = hexagons_yr %>%
 left_join(hex_freq_table_yr, join_by(GRID_ID, year)) %>%
 mutate(hom_ct = ifelse(is.na(hom_ct), 0, hom_ct))

Creating a list of hexagon maps:

hex_map_list = list()

for (n in 1940:1964) {
  data = hexagons_hom_ct_yr %>% filter(year == n)
  map = tm_shape(data) + tm_polygons(col = "hom_ct", breaks = c(0, 1, 5, 10, Inf), title = paste("Homicides in", n))
  hex_map_list[[n-1939]] = map
  rm(n, data, map)
}

Finally, the tmap_animation function is used to create the animation presented below:

tmap_animation(hex_map_list, "hex_map_list.gif", width = 400, height = 500, delay = 50)

Animation of Homicides per Hexagon by Year

Community Areas

A similar set of code is used to create an animation of homicides per Community Area for each year:

# Create table of frequency by CA and year
CA_freq_table_yr = as.data.frame(table(homicides_CA_joined$area_num_1, homicides_CA_joined$year)) %>%
  rename(area_num_1 = Var1, year = Var2, hom_ct = Freq)

# Convert factors to integers/strings
CA_freq_table_yr$area_num_1 = as.character(paste(CA_freq_table_yr$area_num_1))
CA_freq_table_yr$year = as.integer(paste(CA_freq_table_yr$year))

# Create an expanded grid of CAs that includes each CA/year combo
CAs_yr = expand.grid(area_num_1 = CAs$area_num_1, year = 1940:1964) %>%
  left_join(CAs, by = "area_num_1") %>%
  st_as_sf()

# Merge frequency table back into CA/year grid
CAs_hom_ct_yr = CAs_yr %>%
  left_join(CA_freq_table_yr, join_by(area_num_1, year)) %>%
  mutate(hom_ct = ifelse(is.na(hom_ct), 0, hom_ct))

Creating a list of community area maps:

CA_map_list = list()

for (n in 1940:1964) {
  data = CAs_hom_ct_yr %>% filter(year == n)
  map = tm_shape(data) + tm_polygons(col = "hom_ct", breaks = c(0, 1, 5, 10, 50, Inf), title = paste("Homicides in", n))
  CA_map_list[[n-1939]] = map
  rm(n, data, map)
}

Finally, the tmap_animation function is used to create the animation presented below:

tmap_animation(CA_map_list, "CA_map_list.gif", width = 400, height = 500, delay = 50)

Animation of Homicides per Community Area by Year

Preliminary findings: Because homicides are somewhat sporadic events over the course of a 1-year timespan, there appears to be a lot of oscillation between frames. This indicates that a longer timespan should be used for each frame. The animations below present the results with a 5-year timespan instead of 1-year.

5-year Time Series Choropleth Maps

hex_map_list_5yr = list()
CA_map_list_5yr = list()

for (n in seq(1940, 1964, 5)) {
  data = hexagons_hom_ct_yr %>%
    filter(between(year, n, n+4))
  aggregated_data = aggregate(data$hom_ct, list(data$GRID_ID), FUN = sum) %>%
    rename(GRID_ID = Group.1, hom_ct = x) %>%
    left_join(hexagons, by = "GRID_ID") %>%
    st_as_sf()
  map = tm_shape(aggregated_data) + tm_polygons(col = "hom_ct", breaks = c(0, 1, 5, 10, 50, Inf), title = paste0("Homicides in ", n, "-", n+4))
  hex_map_list_5yr[[1+((n-1940)/5)]] = map
  rm(n, data, aggregated_data, map)
}

for (n in seq(1940, 1964, 5)) {
  data = CAs_hom_ct_yr %>%
    filter(between(year, n, n+4))
  aggregated_data = aggregate(data$hom_ct, list(data$area_num_1), FUN = sum) %>%
    rename(area_num_1 = Group.1, hom_ct = x) %>%
    left_join(CAs, by = "area_num_1") %>%
    st_as_sf()
  map = tm_shape(aggregated_data) + tm_polygons(col = "hom_ct", breaks = c(0, 1, 5, 10, 50, 100, Inf), title = paste0("Homicides in ", n, "-", n+4))
  CA_map_list_5yr[[1+((n-1940)/5)]] = map
  rm(n, data, aggregated_data, map)
}

tmap_animation(hex_map_list_5yr, "hex_map_list_5yr.gif", width = 400, height = 500, delay = 100)
tmap_animation(CA_map_list_5yr, "CA_map_list_5yr.gif", width = 400, height = 500, delay = 100)

Preliminary findings: Now the outward spread of homicides is a bit clearer, especially in the CA animation above. However, some oscillation is still present.

Directional Distribution

The following code creates and visualizes the standard deviational ellipse for the homicide data:

sde_data = std_dev_ellipse(homicides)
sde = st_ellipse(geometry = sde_data, sx = sde_data$sx, sy = sde_data$sy, rotation = -sde_data$theta)

leaflet(homicides_jitter) %>%
  addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
  addGlPoints(homicides_jitter, radius = 2, fillColor = "black") %>%
  addPolygons(data = sde, opacity = 0.2)

Preliminary findings: Rather intuitively, the standard deviational ellipse encompasses the South and West Sides, with the Near North Side included as well. Future analysis could demonstrate the change of the ellipse’s size and shape over time, which could demonstrate how the geographic distribution of homicides spread out in a more systematic way.

Part 2: Map Visualizations for Cook County Homicides, 1874-2016

Loading Data

A frequency table of homicides per hexagon before and after the time period shown above, compiled previously, is loaded in.

# Load in homicides per hexagon data
pt2_hex_freq_table_yr = read.csv("homicides_1874-2020_hex/homicides_pre_1931_post_1964.csv") %>%
  select(GRID_ID, year, hom_ct = homicides)

Choropleth Maps

First, all homicides 1874-1930 and 1965-2020 are visualized using a simple choropleth map. Note that records between 1874-1930 include all homicides; records used in the above analysis of homicides between 1940-1964 only include murders, justifiable homicides, and cases in which the verdict is unclear; and records between 1965-2016 only include murders. Thus, the findings from each of these time periods are not directly comparable to each other.

# Separate data between 1874-1930 and between 1965-2020
pt2_before_only = pt2_hex_freq_table_yr %>%
  filter(year < 1931) %>%
  select(GRID_ID, hom_ct)

pt2_after_only = pt2_hex_freq_table_yr %>%
  filter(year > 1964) %>%
  select(GRID_ID, hom_ct)

# Take aggregate of tables by GRID_ID
pt2_before_aggregated = aggregate(.~GRID_ID, data = pt2_before_only,FUN = sum)
pt2_after_aggregated = aggregate(.~GRID_ID, data = pt2_after_only, FUN = sum)

# Merge frequency table into hexes
pt2_before_hexagons_hom_ct = hexagons %>%
  left_join(pt2_before_aggregated, by = "GRID_ID")

pt2_after_hexagons_hom_ct = hexagons %>%
  left_join(pt2_after_aggregated, by = "GRID_ID")

Viewing the hexagons:

tmap_mode("view")
tm_shape(pt2_before_hexagons_hom_ct) + tm_polygons(col = "hom_ct", breaks = c(0,1,5,10,50,100,Inf), title="Homicides 1874-1930", id="hom_ct", popup.vars = c("Grid ID:" = "GRID_ID", "Homicides:" = "hom_ct"))

tm_shape(pt2_after_hexagons_hom_ct) + tm_polygons(col = "hom_ct", breaks = c(0,1,5,10,50,100,Inf), title="Homicides 1965-2020", id="hom_ct", popup.vars = c("Grid ID:" = "GRID_ID", "Homicides:" = "hom_ct"))

The same custom classification scheme described above is used here for consistency, though without the “500 or more” bin because it is not used.

Preliminary findings: These two maps further demonstrate how homicides have spread outward over time, well beyond the time period of 1940-1964.

Time Series Choropleth Map

Next, time series animations like the ones shown above are created for hexagons between 1874-1930 and 1965-2016. A 5-year interval is used like before to reduce noise. Because of this, the years 1874, 1930, 2015, and 2016 will not be included, as they do not fit into any 5-year intervals.

# Create an expanded grid of hexagons that includes each hex/year combo
pt2_hexagons_yr = expand.grid(GRID_ID = hexagons$GRID_ID, year = c(1874:1930, 1965:2016)) %>%
  left_join(hexagons, by = "GRID_ID") %>%
  st_as_sf()

# Merge into hexagons
pt2_hexagons_hom_ct_yr = pt2_hexagons_yr %>%
  left_join(pt2_hex_freq_table_yr, join_by(GRID_ID, year))

Creating a list of hexagon maps:

pt2_map_list_before = list()
pt2_map_list_after = list()

for (n in seq(1875, 1929, 5)) {
  data = pt2_hexagons_hom_ct_yr %>%
    filter(between(year, n, n+4))
  aggregated_data = aggregate(data$hom_ct, list(data$GRID_ID), FUN = sum) %>%
    rename(GRID_ID = Group.1, hom_ct = x) %>%
    left_join(hexagons, by = "GRID_ID") %>%
    st_as_sf()
  map = tm_shape(aggregated_data) + tm_polygons(col = "hom_ct", breaks = c(0, 1, 5, 10, 50, Inf), title = paste0("Homicides in ", n, "-", n+4))
  pt2_map_list_before[[1+((n-1875)/5)]] = map
  rm(n, data, aggregated_data, map)
}

for (n in seq(1965, 2014, 5)) {
  data = pt2_hexagons_hom_ct_yr %>%
    filter(between(year, n, n+4))
  aggregated_data = aggregate(data$hom_ct, list(data$GRID_ID), FUN = sum) %>%
    rename(GRID_ID = Group.1, hom_ct = x) %>%
    left_join(hexagons, by = "GRID_ID") %>%
    st_as_sf()
  map = tm_shape(aggregated_data) + tm_polygons(col = "hom_ct", breaks = c(0, 1, 5, 10, 50, Inf), title = paste0("Homicides in ", n, "-", n+4))
  pt2_map_list_after[[1+((n-1965)/5)]] = map
  rm(n, data, aggregated_data, map)
}

Finally, the tmap_animation function is used to create the two animations presented below:

tmap_animation(pt2_map_list_before, "hex_map_list_before.gif", width = 400, height = 500, delay = 100)
tmap_animation(pt2_map_list_after, "hex_map_list_after.gif", width = 400, height = 500, delay = 100)

Preliminary findings: The two animations above further demonstrate the spreading-out of homicides over time that parallels the outward growth of the city. The first animation, including homicides from 1865 to 1930, begins with homicides almost entirely concentrated in the city center. By 1925-1929, however, homicides are just as concentrated in the Near West and Near South Sides as they are in the city center. The latter animation also shows a diffusion of homicides over time, though also shows a decrease in their overall concentration within the last two time periods. It also suggests that the shape of this distribution has fundamentally changed over the past few decades – from a hub-and-spoke shape centered on the Loop and reaching into the South, West, and North Sides in the 1960s to two distinct regions on the West and South Sides and no central concentration in the Loop in the 2010s.

Time Series Choropleth Map - 1874-2014

An animation of homicides per hexagon per 5-year interval for all time periods is created below.

tmap_animation(c(pt2_map_list_before, hex_map_list_5yr, pt2_map_list_after), "hex_map_list_all.gif", width = 400, height = 500, delay = 100)

Animation of Homicides per Hexagon by 5-Year Timespan 1874-2014

Note: As mentioned above, methodologies varied between the three time periods (all homicides were included in 1874-1930; only murders, justifiable homicides, and unclear verdicts in 1940-1964; and only murders in 1965-2016) and thus time periods should be compared with caution.

Part 3: Weighting Yearly Homicide Counts by Building Counts

In this section, the map of homicide counts per hexagon is reproduced, but with Zillow “ZTRAX” building counts used as a denominator, to get a rate of homicides based on how settled the area was at the time. An additional time series map animation of building counts is provided as well for reference.

hex_bldgs = read.csv("hexagons/hex_building_counts.csv")

pt2_hexagons_hom_ct_yr_ALL = rbind(hexagons_hom_ct_yr, pt2_hexagons_hom_ct_yr)

pt2_map_list_bldgs = list()
pt2_map_list_bldg_weighted = list()

for (n in seq(1875, 2014, 5)) {
  if (n < 1930 | n > 1939) {
    # Homicide data
    data = pt2_hexagons_hom_ct_yr_ALL %>%
      filter(between(year, n, n+4))
    aggregated_data = aggregate(data$hom_ct, list(data$GRID_ID), FUN = sum) %>%
      rename(GRID_ID = Group.1, hom_ct = x) %>%
      left_join(hexagons, by = "GRID_ID") %>%
      st_as_sf()

    # Building data
    bldgs = hex_bldgs %>%
      filter(between(year, n, n+4))
    aggregated_bldgs = aggregate(bldgs$weighted_building_count, list(bldgs$GRID_ID), FUN = mean) %>%
      rename(GRID_ID = Group.1, bldg_ct = x)

    # Combined data
    aggregated_data = aggregated_data %>%
      left_join(aggregated_bldgs, by = "GRID_ID") %>%
      mutate(hom_rt = hom_ct / bldg_ct) %>%
      mutate(hom_rt = ifelse(hom_rt == Inf, NA, hom_rt))

    # Maps
    bldgs_map = tm_shape(aggregated_data) + tm_polygons(col = "bldg_ct", breaks = c(0,1,5,10,50,100,500,1000,Inf), title = paste0("Buildings in ", n, "-", n+4))
    map = tm_shape(aggregated_data) + tm_polygons(col = "hom_rt", breaks = c(0, 0.001, 0.01, 0.1, 1, 10), title = paste0("Homicides per building in ", n, "-", n+4))

    pt2_map_list_bldgs[[1+length(pt2_map_list_bldgs)]] = bldgs_map
    pt2_map_list_bldg_weighted[[1+length(pt2_map_list_bldg_weighted)]] = map
    
    rm(n, data, aggregated_data, bldgs, aggregated_bldgs, map, bldgs_map)
  }
}

An animation of building-weighted homicide rates per hexagon per 5-year interval is presented below:

tmap_animation(pt2_map_list_bldg_weighted, "hex_map_list_bldg_weighted.gif", width = 400, height = 500, delay = 100)

Animation of Homicides per Building per Hexagon by 5-Year Timespan 1874-2014

Note: Areas labeled “missing” had no buildings, so the calculation of the rate produced a division by 0 error.

Click here to see a time series animation of building counts

tmap_animation(pt2_map_list_bldgs, "hex_map_list_bldgs.gif", width = 400, height = 500, delay = 100)

Animation of Buildings per Hexagon by 5-Year Timespan 1874-2014

Part 4: Time Series Line Graphs by Annexation Status and Land Sale Year

In this section, time series line graphs are made to show the rate of homicides per building for each year in hexagons categorized by annexation status and land sale year.

Creating a Time Series Dataset with Homicides per Building Rates for All Hexagons and Years

Previous homicide and building data are used to make a dataset for Part 3 that characterizes each hexagon by its homicide count, building count, and homicide per building rate for each year.

pt4_data = left_join(pt2_hexagons_hom_ct_yr_ALL, hex_bldgs, by = c("GRID_ID", "year")) %>%
  select(GRID_ID, year, hom_ct, bldg_ct = weighted_building_count) %>%
  mutate(hom_rt = hom_ct / bldg_ct) %>%
  mutate(hom_rt = ifelse(hom_rt == Inf, NA, hom_rt))

Annexations Status

First, hexagons are categorized by the following classification scheme: * Hexagons within the original Chicago city limits * Hexagons within areas incorporated into the city by the state * Hexagons within settlements that were annexed into the city * Hexagons within settlements (Jefferson and Lake Townships) that were sparsely populated when they were annexed into the city.

This is done by calculating the centroids of the hexagons, and spatially assigning those centroids to the surrounding area. Additionally, any hexagon centroids that are outside of the city limits are excluded.

annexed_areas = read_sf("annexed/annexed_areas/Annexed Settlements.shp") %>%
  st_transform(st_crs(hexagons))

hexagons_annexed = hexagons %>%
  st_centroid() %>%
  st_filter(city_bounds, .predicate = st_within) %>%
  st_join(annexed_areas, join = st_within) %>%
  select(-year) %>%
  mutate(annex_category = case_when(
    settlement == "CHICAGO" ~ "Original",
    is.na(settlement) ~ "Incorporated",
    settlement == "LAKE" | settlement == "JEFFERSON" ~ "Annexed - Sparse",
    .default = "Annexed - Populated"
  )) %>%
  st_drop_geometry() %>%
  right_join(hexagons, ., by="GRID_ID") %>%
  select(GRID_ID, settlement, annex_category)

The map below summarizes the categorization of hexagons by annexation status:

tmap_mode("view")
tm_shape(hexagons_annexed) + tm_polygons(col = "annex_category", title="Annexation Category", id="annex_category", popup.vars = c("Grid ID:" = "GRID_ID", "Settlement:" = "settlement", "Annexation Category:" = "annex_category"))

Next, the annexation table is joined to the Part 3 dataset to categorize each observation by annexation status. Then, observations are aggregated by annexation status to produce the time series graph shown below.

pt4_annexed = right_join(pt4_data, st_drop_geometry(hexagons_annexed), by = "GRID_ID") %>%
  st_drop_geometry()

pt4_annexed_aggregated = pt4_annexed %>%
  group_by(year, annex_category) %>%
  summarise(hom_ct = sum(hom_ct), bldg_ct = sum(bldg_ct)) %>%
  mutate(hom_rt = hom_ct / bldg_ct)

# Line Graph
ggplot(data = pt4_annexed_aggregated, aes(x = year, y = hom_rt, group = annex_category)) +
  geom_line(aes(color = annex_category), size = 1) +
  ylab("Homicides per Building") +
  xlab("Year") +
  coord_cartesian(ylim = c(0, 0.01))

Note: The y-axis is limited at 0.01 to reduce the impact of noise on the interpretability of the rest of the graph.

Map Animation: Homicide Rate by Year by Annexation Category

The animation below shows homicide rates by year, separated by annexation category.

mean_na_rm = function(x) {
  mean(x, na.rm = TRUE)
}

pt4_map_list_annexed = list()

for (n in seq(1875, 2014, 5)) {
  if (n < 1930 | n > 1939) {
    data = pt4_annexed %>%
      filter(between(year, n, n+4))
    aggregated_data = aggregate(data$hom_rt, list(data$GRID_ID), FUN = mean_na_rm) %>%
      rename(GRID_ID = Group.1, hom_rt = x) %>%
      left_join(hexagons_annexed, by = "GRID_ID") %>%
      st_as_sf()

    map = tm_shape(aggregated_data) + tm_facets(by = "annex_category", free.coords = F, scale.factor = 2) + tm_polygons(col = "hom_rt", breaks = c(0, 0.001, 0.01, 0.1, 1, 10), title = paste0(n, "-", n+4))

    pt4_map_list_annexed[[1+length(pt4_map_list_annexed)]] = map
    
    rm(n, data, aggregated_data, map)
  }
}

An animation of building-weighted homicide rates per hexagon per 5-year interval is presented below:

tmap_animation(pt4_map_list_annexed, "hex_map_list_annexed.gif", width = 800, height = 1000, delay = 100)

Animation of Homicides per Building per Hexagon by 5-Year Timespan 1874-2014; Separated by Annexation Category

Land Sale Year

Work in progress

Redlining

Finally, hexagons are categorized by redlining grades, A, B, C, D, and other (includes industrial, commercial, or sparsely built-up areas).

This is done by calculating the centroids of the hexagons, and spatially assigning those centroids to the surrounding area’s redlining grade. Additionally, any hexagon centroids that are outside of the city limits are excluded.

redlined_areas = read_sf("redlining/redlining.geojson") %>%
  st_transform(st_crs(hexagons)) %>%
  st_make_valid() %>%
  select(grade)

hexagons_redlined = hexagons %>%
  st_centroid() %>%
  st_filter(city_bounds, .predicate = st_within) %>%
  st_join(redlined_areas, join = st_within) %>%
  st_drop_geometry() %>%
  right_join(hexagons, ., by="GRID_ID") %>%
  select(GRID_ID, grade)

The map below summarizes the categorization of hexagons by annexation status:

tmap_mode("view")
tm_shape(hexagons_redlined) + tm_polygons(col = "grade", palette = c("#00aa00", "#0000ff", "#ffff00", "#ff0000"), title="Redlining Grade", id="grade", popup.vars = c("Grid ID:" = "GRID_ID", "Redlining Grade:" = "grade"))

Next, the redlining table is joined to the Part 3 dataset to categorize each observation by redlining grade. Then, observations are aggregated by redlining grade to produce the time series graph shown below.

pt4_redlined = right_join(pt4_data, st_drop_geometry(hexagons_redlined), by = "GRID_ID") %>%
  st_drop_geometry()

pt4_redlined_aggregated = pt4_redlined %>%
  group_by(year, grade) %>%
  summarise(hom_ct = sum(hom_ct), bldg_ct = sum(bldg_ct)) %>%
  mutate(hom_rt = hom_ct / bldg_ct)

# Line Graph
ggplot(data = pt4_redlined_aggregated, aes(x = year, y = hom_rt, group = grade)) +
  geom_line(aes(color = grade), size = 1, palette = c("#00aa00", "#0000ff", "#ffff00", "#ff0000")) +
  ylab("Homicides per Building") +
  xlab("Year") +
  scale_color_manual(values = c("#00aa00", "#0000ff", "#ffff00", "#ff0000")) +
  coord_cartesian(ylim = c(0, 0.01))

Note: The y-axis is limited at 0.01 to reduce the impact of noise on the interpretability of the rest of the graph.

Map Animation: Homicide Rate by Year by Annexation Category

The animation below shows homicide rates by year, separated by annexation category.

pt4_map_list_redlined = list()

for (n in seq(1875, 2014, 5)) {
  if (n < 1930 | n > 1939) {
    data = pt4_redlined %>%
      filter(between(year, n, n+4))
    aggregated_data = aggregate(data$hom_rt, list(data$GRID_ID), FUN = mean_na_rm) %>%
      rename(GRID_ID = Group.1, hom_rt = x) %>%
      left_join(hexagons_redlined, by = "GRID_ID") %>%
      st_as_sf()

    map = tm_shape(aggregated_data) + tm_facets(by = "grade", free.coords = F, scale.factor = 2) + tm_polygons(col = "hom_rt", breaks = c(0, 0.001, 0.01, 0.1, 1, 10), title = paste0(n, "-", n+4))

    pt4_map_list_redlined[[1+length(pt4_map_list_redlined)]] = map
    
    rm(n, data, aggregated_data, map)
  }
}

An animation of building-weighted homicide rates per hexagon per 5-year interval is presented below:

tmap_animation(pt4_map_list_redlined, "hex_map_list_redlined.gif", width = 800, height = 1500, delay = 100)

Animation of Homicides per Building per Hexagon by 5-Year Timespan 1874-2014; Separated by Redlining Grade

Part 5: Survivorship Graph of Hexagons with Few Homicides per Year

Work in progress